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ABSTRACT 

We study the impact of ultraviolet background (UVB) radiation field and the local stellar radiation 
on the Hi column density distribution f{Nm) of damped Lya systems (DLAs) and sub-DLAs at 
z = 3 using cosmological smoothed particle hydrodynamics simulations. We find that, in the previous 
simulations with an optically thin approximation, the UVB was sinking into the H i cloud too deeply, 
and therefore we underestimated the /(A^hi) at 19 < logA'ni < 21.2 compared to the observations. 
If the UVB is shut off in the high-density regions with rigas > 6 x 10~^cm~^, then we reproduce 
the observed /(iVni) at 2; = 3 very well. We also investigate the effect of local stellar radiation by 
post-processing our simulation with a radiative transfer code, and find that the local stellar radiation 
does not change the /(iVni) very much. Our results show that the shape of /(A^hi) is determined 
primarily by the UVB with a much weaker effect by the local stellar radiation and that the optically 
thin approximation often used in cosmological simulation is inadequate to properly treat the ionization 
structure of neutral gas in and out of DLAs. Our result also indicates that the DLA gas is closely 
related to the transition region from optically-thick neutral gas to optically-thin ionized gas within 
dark matter halos. 

Subject headings: cosmology: theory — galaxies: evolution — galaxies: formation — galaxies: high- 
redshift — quasars: absorption lines — methods: numerical 



1. INTRODUCTION 

The Hi column density distribution function /{Nm) 
is one of the most basic statistics of quasar absorption 
systems, similarly to the luminosity function of galaxies. 
The accuracy of observational data on /(A'hi) has dra- 
matically improved over the past several years, thanks to 
the large samples of damped Lya systems (DLAs) and 
sub-DLA s discovered in large data sets of quasar spec- 
tra (e.g., Pcroux et al.' 2003; 'Prochaska & Herbert-FortI 
2004 : Pcroux et al. .2005; Prochaska ct al. 2005, 2008; 



Noterdaeme et al.ll2009l) . Their results indicate that the 
observed /(A^hi) of DLAs can be fitted well with either 
a double power-law or a Schechter-type function. 

It would be desirable to understand the physical ori- 
gin of the shape of /{Nhi) in a cosmological con- 
text of the standard A cold dark matter mo d el. Th e 
first such attempt was made by iKatz et al.l (|1996b( ). 
who used a cosmological smoothed particle hydrody- 
namics (SPH) simulation with a comoving box-size of 
22.22/1^^ Mpc, 2 X 64"^ particles, and cosmological param- 
eters {nm,^A,^b) = (1.0,0.0,0.05). They found that 
their simulation underpredicted /(A^hi) compared to the 
observational data by a factor of few or more, using a 
uniform ultraviolet background (UVB) J(z/) = Jo(fo/i^) 
at z = 3, where z/q is the Lyman-limit frequency and 
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Jo = 10-22 ergs-1 cm-2 sr'^ Hz"!. The UVB is usually 
treated with an optically thin limit regardless of gas den- 
sity in cosmological hydrodynamic simulations owing to 
the computational limit. This simplified approximation 
may artificially increase the ionization fraction of gas, 
and gives r ise to t he discrepancy in /(A^hi)- 

.Nagamine et al.l (2004, 2007) updated the work of 
IKatz et all (|1996bD using cosmological SPH simulations 
with a comoving box size of 10h~^ Mpc, 2x324^ particles, 
and {nm,^A,^b) = (0.3,0.7,0.044). Interestingly, they 
also found a similar underprediction of /(A^hi) compared 
to the observatio ns, despite signific antly higher resolu- 
tion than that of iKatz et al.l (Il996b[). Their simu l ations 
included a uniform U VB of iHaardt fc Madaul (|1996D 
spectrum, modified by iDave et al.l (|1999( ) to match the 
Lya forest observations. 

We have attempted to resolve this discrepancy by mod- 
ifying the models of star formation (SF) and supernova 
feedback; e.g., changing the SF threshold density, SF 
tim e-scale, feedback strengths , or adding metal-line cool- 
ing (|Choi &: Nagamin3l2009bl laf) . However, none of these 
changes in the physical models resolved the discrepancy 
in /(A'hi) fundamentally. 

In this Letter, we show that the effect of UVB is the key 
in determining the shape of /(A'hi) at logA^Hi ^ 21.6. 
In addition, we consider the effect of local stellar radia- 
tion on /(A^hi) by performing a radiative transfer (RT) 
calculation. Our paper is organized as follows. In Sec- 
tion [21 we briefly describe the setup of our simulations, 
and present the results in Section [3l We then discuss 
the comparison with other recent works, and conclude in 
Section 111 

2. SIMULATIONS 

We use the updated version of the tree-parti cle- mesh 
SPH code GADGET-3 (originally described in iSpringell 
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120051) . Our convention al code inc ludes ra diative cooling 
by H, He, and metals ()Choi fc N agamin d l2009aD, heat- 
ing by a uniform UVB of a modified H aardt fc Madaul 
(| 19961 ) spectrum (|Katz et al.l [T996a: Da ve et al.l 119991 ) . 
SF, supernova feedback, a phenomenological model for 
galactic winds, and a sub-r esolution model of multiphas e 
interstellar medium (ISM; iSpringel fc Hernquisti 120031 ) . 
In this multiphase ISM model, high-density ISM is pic- 
tured to be a two-phase fluid consisting of cold clouds 
in pressure equilibrium with a hot ambient phase. Cold 
clouds grow by radiative cooling out of the hot medium, 
and this material forms the reservoir of baryons avail- 
ab le for SF. We use the "P ressure SF" model described 
by Choi fc Nagamind (120101 ) (whi ch is based on the work 
bv lSchave fc Dalla Vecchiar()2008[ )). but we have checked 



that the details of the SF model do not change the main 
conclusions of this paper. 

For all the simulations used in this paper, we employ 
a box size of comoving 10 h"^ Mpc and a total particle 
number of 2 x 144^ for gas and dark matter. The initial 
gas particle mass is mgas = 4.1 x 10^ h~^MQ, and the 
dark matter particle mass is mdm = 2.0 x 10^ h~^MQ. 
The comoving gravitational softening length is 
2.78/i^^kpc, so the physical resoluti on of our sim- 
ulatio n is ~0.7/i~^kpc at z = 3. iNagamine et al.l 
(|2004D showed that increasing the particle number from 
2 X 144^ to 2 x 324^ did not change the shape of /(Nm) 
very much, therefore our results would not be strongly 
affected by the resolution effect. (But see further discus- 
sion in SectionlD) The comoving box size of 10 Mpc 
is somewhat small, however, the number of missed very 
massive haloes are relatively small, and the imp act on 
/£/Vm) is expected to be small. In fact , Nagami ne et al.l 
' (|2004l ) showed that /(Nm) did not change very much 
with increasing box size, except that the lower iVui end 
of fi^m) decreased due to lower resolution. In this 
work, we have not corrected our results for the box 
size effect. The adopted cosmological parameters of all 
simulations are consis t ent w ith the latest WMAP result 
(|Komatsu et al.l 120091 l2010( ): (n^, Qa, ^h, cs, h, n^) = 
(0.26,0.74,0.044,0.80,0.72,0.96), where h 
iJo/(100kms-iMpc-i). 

With this setup, we run four simulations with different 
models of UVB: "Fiducial", "No-UV", "Half-UV", and 
"OTUV" (Optically Thick UV) runs. In the Fiducial 
run, the gas is heated and ionized by the uniform UVB 
under the optically thin approximation. In the No-UV 
run, the UVB strength is set to zero. In the Half-UV 
run, the normalization of UVB is reduced by half. In the 
OTUV run, we assume that the uniform UVB cannot 
penetrate into the high-density gas with rigas > n 
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but otherwise it is the same as the Fiducial run at Ugas < 
„uv 

We adopt the threshold density n 
6 X 10~^cm~^, where n^jf is the SF threshold density 
above which the stars are allowed to form. In our simu- 
lations, the gas with rigas > nf^ is mostly neutral owing 
to the multiphase ISM model. We originally arrived at 
the above value of n^-^^ by successively lowering its value 
from n^^ and checking the agreement with the observed 
f{Niii), but will provide further justifications below. 

There is evidence that the above value of n^^^ is phys- 
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Fig. 1. — Hi column density distribution functions at 2 = 3 for 
the four runs with different treatmen t of U VB. The observational 
data points are from Pcroux et al? f2005', black open squaresV 
O'Meara et aLj I I2007I . magenta open circles), Prochaska & Wol^ 
(|2009i . green triangles) , and .Noterdaeme et al.. (,200S . cyan bars) . 



ically appropriate. iTaiiri fc Umemural ()1998[ ) found that 
the hydrogen cloud becomes fully self-shielded above a 
critical density of 1.4 x lO^'^ cni^'^ through RT calcula- 
tions for a spherical top-hat sphere, and that the critical 
density has a mild depende nce on t he clo ud mass and 
the UVB intensity. KoUmci er et al.l (|2010f ) performed a 
three-dimensional UVB RT calculation with an isother- 
mal sphere, and showed that the above value of a,p- 
proximatel y corresponds to the transitio n density from 
Hii to Hi. iFaucher-Giguere et al.l (|20100 postprocessed 
cosmological SPH simulations with a ray tracing code 
and found that turning off UVB at rigas > 0.01 cm~'^ 
produces a favorable result. Furthermore, we also con- 
firmed that the above is appropriate by postprocess- 
ing our simulations with a RT code, which we will report 
in detail in a separate paper (H. Yajima et al., 2011, 
in preparation). For these reasons, we consider that the 
correct value ofn^^ is in the range of 10~^ to 10^"^ cm^'^, 
depending on the cloud mass and UVB intensity. 

3. RESULTS 

3.1. Effect of UVB on /(iVm) 

Figure[T] shows the /(A^m) in the four runs with differ- 
ent UVB treatment, which was calcu l ated b y the same 
method described in INagamine et alj ()2004[ ). In short, 
we set up a uniform grid around each dark matter halo, 
and project the gas density field onto a face of the grid to 
compute A^Hi- Figure[T] clearly shows that the Fiducial 
run underpredicts the /(Ahi), particularly at logA^Hi < 
21.2. The Half-UV run is somewhat higher than the 
Fiducial run, but still not enough to account for the ob- 
served number of columns at 19.5 < logiVni < 21. On 
the other hand, the No-UV run completely overpredicts 
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Fig. 2. — Neutral hydrogen fraction (ngas vs. Xhi; top panels) and phase space distribution (rigas vs. T; bottom panels) of the cosmic 
gas at 2 = 3 for the four runs. The neutral fraction is defined as xhi = ^Hl/('^HI + "Hn)- For plotting purpose, we plot only randomly 
selected 1% of the total gas particles in the simulation. The two densities, 
and blue dotted lines, respectively. 

the observed /(iVni) at all column densities. The sudden 
decrease in f{Nm) from the No-UV run to the Half-UV 
run shows that even a weak UVB can ionize the gas and 
the simulation cannot account for the observed number 
of columns. The OTUV run agrees very well with the 
observed data at 19.2 < logA'ni < 21.2. The compari- 
son of these four runs suggests that the UVB affects the 
shape of /(A^hi) significantly, and that applying the opti- 
cally thin approximation at all densities is too simplistic 
to reproduce the observed /(-/Vhi) properly. Our simu- 
lations overpredict fiNm) at logA^ni > lO^^cm"^. But 
this could be because our simulations do not take into 
account the conversion of H i into H2 , and the conclusion 
of this paper does not depend on this issue. 

To obtain a better physical intuition on the effect of 
UVB, we plot the hydrogen neutral fraction (xhi) and 
temperature against gas density for the four runs in Fig- 
ure [5] The difference of /(A^hi) between the Fiducial 
and the OTUV run must be due to the different neutral 
fraction of gas at < rigas < n^^ in the two runs. 
The Fiducial run predicts a high degree of ionization at 
rigas < n^^ (« 0.6 cm~^). For example, the gas with 
rigas ~ 0.1 cm~'^ at z = 3 has xhi ~ 0.05 for the Fiducial 
run. By referencing to /(A'hi), the gas with this density 
should correspond to the DLAs with logiVni ~ 20, and 
Xhi ~ 0.05 is too low to match the observed data on 
f{Nm)- The Half-UV run has a higher neutral fraction 
for the gases with ngas < n-^^, but still with xhi < 0.3. In 
contrast, most hydrogen in the No-UV run has xhi ~ 1-0, 
except for the small fraction of hot gas which is ionized 
by collisional ionization in shocked regions. Eliminating 
the UVB completely is not a realistic assumption, and 
the No-UV run clearly overpredicts /{Nhi) as expected. 
Lastly, in the OTUV run, the gas with rigas > has 
Xhi ~ 1-0, which allows it to match the observed /(TVhi). 
Figures[T] and [5] clearly show that the optically thin ap- 
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too much. In order to properly reproduce the observed 
f{Nm) with our simulations, the gas should be mostly 
neutral at rigas > n^]^ . 

The effects of UVB can also be seen in the p — T phase 
diagrams in the bottom panels of Figure [51 The com- 
parison of the Fiducial and No-UV runs shows the well- 
known photoionization effect on the diffuse IGM at low 
densities of rigas < 10~^cm~^. Compared to the Fidu- 
cial run, the gas with n]^^ < rigas < "-tif i'^ the OTUV 
run is almost fully neutral with a lower temperature of 
T ^ 10** K. We will show in a separate paper (H. Yajima 
et al., 2011, in preparation) that a full RT calculation 
supports the results shown in FigurelH^d). 

3.2. Effects of Radiative Transfer on /(iVni) 

The OTUV run is remarkably successful in reproducing 
the observed /(-/Vri) at 19 < logTVni < 21.5, but so 
far we have neglected the radiative effects of the stellar 
radiation from local stellar sources, which could heat and 
ionize the high-density gas in the star-forming regions. 

In order to consider the effects of local stellar radia- 
tion, we postprocess our simulati on with the Authentic 
Radiation Transfer (ART) code ()Nakamoto et all 120011: 
lYaiima et al.l I2009D . which uses the ray-tracing tech- 
nique. We have performed t he RT calc ulation in a similar 
fashion to the one we did in lYaiima e t al. (2010), taking 
the RT grid cell size equal to the gravitational softening 
length. This way the resolution of the RT calculation is 
fixed for all halos. Figure [3] compares the results from 
the OTUV run with and without the RT calculation. 
The transferre d stellar spectra are computed by using the 
PEGASE v2.0 (jFioc fc Rocca-Volmerangell997D based on 
the mass, formation time, and age of the star particles 
generated in the simulation. We find that the local stel- 
lar radiation does not have a significant impact on the 
shape of f{NHi), as shown by the red long-dashed line 
in Figure[3l which is almost overlapping with the origi- 
nal OTUV run without RT. This is because the stellar 
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Fig. 3. — Effect of local stellar radiation on /(Afjji) at 2 = 3. 
The observational data points are the same as in Figure [l] The 
two lines for the OTUV with and without RT (red dashed and blue 
solid lines, respectively) are almost overlapping with each other. 



radiation can only increase the ionization fraction of gas 
near the ionization boundary where it makes a transi- 
tion from neutral to highly ionized, and the high density 
neutral gas remains largely intact. We will report more 
details of the RT calculation results in a separate publi- 
cation (H. Yajima et al. 2011, in preparation). 

4. DISCUSSIONS AND CONCLUSIONS 

Using cosmological SPH simulations, we examined the 
effects of UVB on /(A^hi), and clarified the reason why 
earlier simulations have underestimated this quantity 
compared to the observations. We find that the radi- 
ation sinks into the halo gas too deeply under the op- 
tically thin approximation of UVB, and the gas with 
19 < logiVni < 21.5 is overly ionized. If we turn off the 
UVB above a physical density of = 6 x 10"'^ cm~^, 
then we reproduce the observed /(iVHi) at z = 3 very 
well. The exact value of would depend on the de- 
tails of the spectral sha pe and the inte nsity of UVB, and 
the size of gas cloud (Taii ri fc Umem ura 1998). Based 
on the comparison to other works, we consider that a 
value of ~ 10~^ — 10~^ cm~^ would be appropriate 
for current cosmological simulations. Our results clearly 
show that the optically thin approximation was respon- 
sible for the failure in matching the observed /(A^hi) in 
earlier simulations. 



Recently, problems in the optically thin approxi- 
mation have been pointed out in the context of the 
Hen reionization eff ect on the thermal history of IGM 
()McQuinn et all [2009 : Fauchcr-GiEuerc e t al.ll2009D and 
Ly-Q emission (Ya ng et al.ll2006t iFaucher-Giguere et al.l 
I2010D . Our present work is another example of in- 
adequacy of optically thin approximation of UVB for 
modeling the ionization of gas near DLAs. The RT 
effects, as well as the exact shape of UVB, have sig- 
nificant impact on the details of galaxy formation his- 
tory including DLAs (e - K-i i Taiiri & Umemura 199J; 
Zheng fc M iralda-Escudc 2002; Hambrick ct all 120091: 
KoUmeier et al..,2010: Faucher-Giguere et al... 20101) . 

Although cosmological RT simulations are beginning 
to be performed, it is still difficult to do a self-consistent 
ray-tracing RT calculation concurrently with the hydro- 
dynamics due to limited computational speed. Under 
this circumstance, it would still help to have a physically 
plausible working model of self-shielding for cosmological 
simulations. Our result provides a useful proxy for the 
threshold density [nY]^) above which the self-shielding 
effect kicks in. 

We also examined the effects of local stellar radiation 
by postprocessing the simulation with a ray-tracing code. 
The effect is not as strong as the UVB, and it only slightly 
decreases /(A^hi)- We will report the further details of 
this ray-tracing work and its effect on DLA cross sec- 
tion in a separate publication (H. Yajima et al. 2011, in 
preparation) . 

To put it in another way, our results suggest that the 
shape of /(A^hi) and the ionization structure of sub- 
DLAs and DLAs would be great probes of UVB at 
high redshift, and that the DLA gas is closely related 
to the transition region from optically-thin ionized gas 
to optically-thick neutral gas within dark matter halos. 
Further comparisons on the physical properties of DLAs, 
sub-DLAs, and Lyman limit systems between simula- 
tions and observations would give us useful insight on 
the nature of UVB. 

In Section[2]we argued that our results are not strongly 
affected by the numerical resolu tion of our curre nt sim- 
ulations, based on the result s oflNagamine et al] ^0041 . 
However the simulations of iNagamine et al.l ()2004l) did 
not include the self- shielding treatment of the OTUV 
run, so to really test the resolution effect, we would have 
to run a new OTUV simulation with 2 x 324'^ particles, 
which we plan to perform in the near future and address 
the resolution effects more directly. 
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